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ABSTRACT 


Progress  is  described  on  the  development  of  an  atmospheric 
fallout  model  to  be  used  to  estimate  particulate  concentrations  in 
the  air  and  on  the  earth's  surface  for  particles  released  at  any  alti¬ 
tude  up  to  90  km.  The  particle  size  range  of  interest  for  this  study 
is  1  to  200  microns  in  diameter.  The  model  is  being  programed  to 
be  run  on  a  CDC  3  600  computer.  During  this  quarter,  preliminary 
specifications  and  coding  for  the  program  were  completed  for  all 
but  one  major  component  of  the  model. 
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SUMMARY 


Knowledge  of  the  fate  of  radioactive  fuel  particles  injected  into  the  atmosphere  by  a 
re-entering  SNAP  device  is  necessary  in  order  to  determine  if  any  hazard  could  exist  to 
the  general  population.  It  is  the  purpose  of  the  work  under  this  contract  to  develop  a  theo¬ 
retical  atmospheric  model,  programmed  to  be  run  on  a  CDC  3600  computer,  that  can  be 
used  to  predict  the  fallout  pattern  and  ultimate  ground  concentrations  of  particulates  in  the 
size  range  of  1  to  200  microns  released  as  a  line  source  anywhere  in  the  atmosphere  up  to 
altitudes  of  90  km. 

In  this  report,  data  is  given  on  the  preliminary  computer  program  specifications  de¬ 
veloped  for  the  fallout  model.  Details  on  the  development  and  coding  of  computer  program 
specifications  for  such  items  as  particle  fall  speeds,  mass-element  trajectories,  as  well 
as  the  methods  for  representing  zonal  flow  and  long  atmospheric  waves  are  presented. 
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PARTICLE  FALLOUT  AND  DISPERSION  IN  THE  ATMOSPHERE 


1.0  INTRODU  CTION 

In  the  I'ir.sl  Quarterly  Progress  Report,  issued  under  Sandia  Corporation 
Contract  No.  4S-2417,  qualitative  descriptions  were  given  of  the  principal  components 
or  subroutines  of  a  computer  program  for  the  fallout  and  dispersion  of  solid  particles 
released  from  a  source  in  the  stratosphere  [8].  The  mathematical  framework  for  the 
simulation  of  large-scale  atmospheric  flow  in  the  model  was  also  presented. 

The  present  report  describes  progress  during  the  second  quarter  (September  — 
November,  1965).  Preliminary  computer  program  specifications  were  completed  for 
the  fallout  model  and  coding  was  completed  for  all  major  components  except  the  mass- 
element  expansion  and  subdivision  routine.  Special  attention  was  given  to  the  problem 
of  selecting  appropriate  particle  fall-speed  equations  for  the  wide  range  of  Reynolds 
numbers,  Knudsen  numbers,  and  Mach  numbers  that  must  be  considered  for  particles 
from  ten  to  sc\'eral  hundred  microns  in  size  falling  from  altitudes  up  to  90  Ion.  The 
transient  acceleration  phase  of  fall  was  investigated  for  particles  with  both  small  and 
large  initial  tall  speeds. 

Two  methods  of  numerical  integration  of  particle  trajectories  were  programmed 
and  tested  using  simulated  equations  of  motion  which  possess  some  of  the  character¬ 
istics  expected  in  more  realistic  trajectories,  and  yet  are  of  sufficiently  elementary 
forms  so  that  exact  solutions  could  be  obtained  for  test  purposes.  The  trajectory  com¬ 
putations  are  one  of  the  most  critical  parts  of  the  model  and  further  testing  is  planned 
for  the  third  quarter. 

Several  computer  programs  required  for  fitting  the  mathematical  functions 
described  in  (S)  to  observed  or  analyzed  winds  were  prepared  and  checked.  Fitting  of 
meridional  profiles  of  the  east— west  or  zonal  mean  winds  at  ten  levels  (summer  and 
winter  seasons)  and  the  large-scale  eddy  flow  at  two  levels  (winter  season)  was  com¬ 
pleted.  Vertical  fitting  of  the  mean  zonal  flow  was  near  completion  at  the  end  of  the 
quarter.  Available  sources  of  data  that  can  be  used  readily  for  this  purpose  have 
apparently  been  exhausted  and  further  progress  will  require  specification  based  on 
known  properties  of  the  flow  and  additional  harmonic  analyses  of  wind  or  pressure 
observations. 
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2.0  PARTICLE  FALLOUT  AND  DISPERSION  MODEL  (Model  B) 


2.1  Source  Description 

Because  Model  B  does  not  take  into  account  the  processes  of  particle  formation, 
it  is  necessary  to  select  one  instant  in  time  during  reentry  that  can  be  used  to  specify 
the  required  initial  conditions  for  the  model.  In  general,  the  elapsed  time  between 
particle  formation  and  the  moment  for  which  initial  conditions  are  specified  should  be 
small  enough  so  that  cloud  growth  by  gravitational  sorting  and  by  turbulent  dispersion 
is  small,  but  large  enough  so  that  their  velocities  are  essentially  independent  of  paiticle 
velocities  at  the  time  of  formation.  The  latter  requirement  is  examined  in  some  detail 
in  Section  2.2.4. 

At  the  initial  moment  (t  =  0  in  Model  B)  it  is  assumed  that  the  particle  cloud  has 
the  form  of  a  simple  line  or  arc  segment  with  a  circular  cross-section  and  finite  length 
roughly  analagous  to  a  meteor  trail  or  a  high-altitude  aircraft  contrail.  In  the  model, 
this  source  cloud  is  simulated  by  one  to  twenty  discrete  elements.  Each  element  takes 
the  form  of  a  trivariate  normal  distribution  of  particles  in  space  characterized  by  the 
following  parameters: 

(a)  Position  (altitude,  latitude,  and  longitude) 

(b)  Dimensions  (standard  deviations  of  the  particle  distribution  in 

the  meridional,  zonal,  and  vertical  directions) 

(c)  Source  Strength  (total  mass  of  particles) 

(d)  Mass-Diameter  Distribution 

In  this  report,  reference  will  be  made  only  to  a  single  source  element,  because  the 
final  concentration  field  and  deposit  patterns  for  the  total  particle  cloud  are  obtained 
by  simple  superposition  and  smoothing  of  contributions  from  all  source  elements.  In 
other  words,  the  transport  and  growth  of  an  individual  source  cloud  element  is  treated 
independently  of  all  other  elements. 

As  indicated  by  (d)  above,  each  source  cloud  element  is  characterized  by  a 
particle-size  distribution.  In  the  model  the  particle  size  distribution  is  divided  into 
discrete  mass  elements,  each  containing  a  fixed  equal  fraction  of  the  total  mass  in  the 
form  of  particles  of  uniform  size.  Three  options  are  available  for  specifying  the  mass- 
diameter  distribution;  a  log-normal  distribution  characterized  by  mass-median 
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'  diameter  6  and  logarithmic  standard  deviation  v;  a  discrete  distribution  character- 

m 

ized  by  mass-median  diameters  for  each  mass  element,  and;  an  unspecified  distribution 
function  other  than  the  log-normal  that  may  be  added  to  the  program  if  desired.  Each 
source  clement  can  be  subdivided  into  1,  20,  50,  or  100  mass  elements.  The  first  of  these 
corresponds  to  a  monodisperse  particle  population.  The  source  description  subroutine 
of  the  computer  program  for  Model  B  has  been  coded  and  verified. 

2.2  Particle  Fall  Speeds 

2.2.1  Atmospheric  Parameters 

In  Model  A,  the  terminal  velocities  of  the  particles  contained  in  each  mass  ele¬ 
ment  were  required  only  at  preselected  altitudes;  consequently,  the  numerical  values 
of  air  density,  viscosity,  and  mean  free  path  were  required  only  at  these  altitudes. 

In  Model  B,  however,  the  altitudes  at  which  these  variables  are  required  are  deter¬ 
mined  during  the  trajectory  integrations.  In  general,  these  altitudes  will  be  unequally 
spaced  and  will  vary  from  one  mass  element  to  another.  For  purposes  of  interpola¬ 
tion,  the  following  poljmomials  were  obtained  by  least-squares  fitting  procedures  using 

-3 

1962  U.S.  standard  atmosphere  values  of  air  density  [p  (gm  cm  )],  viscosity 
[(|j,(gm  cm  ^  sec  ^)],  and  mean  free  path  [X(cm)]  tabulated  at  1.0  km  height  intervals 
from  0.5  km  to  89.5  km  [4]. 

-2  -3  2 

p  =  exp  {  -  6.697  -  8.992  x  10  z  -  2.510  X  10  z 

+  1.206  X  lO^^z^  +  7.748  Xio'^z'^  -  8.199  X  10"®z^ 

”116  ”13T 

-  2.519  X  10  z  +  3.703  X  10  z  }  (2-1) 

”4  —6  — Y  2 

p,  =  1.802  X  10  -  4.362  x  10  z  +  1.173  x  10  z 

-10  3  -11  4 

+  8.877  x  10  z  -  4.715  x  10  z 

-13  5 

+  3.107  x  10  z  (2-2) 

X  =  8.125  x  lO”®  p“^.  (2-3) 
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The  root-mean- square  (rms)  error,  the  mean  absolute  error,  and  the  extreme 
errors  were  all  considered  in  the  selection  of  these  polynomials.  These  quantities  are 
shown  in  Table  2-1  for  polynomials  of  various  degrees. 

TABLE  2-1 

ERROR  SUMMARY  FOR  POLYNOMIAL  REPRESENTATION  OF 
STANDARD  ATMOSPHERE  VALUES  OF  DENSITY  AND 
VISCOSITY  TO  90  km 


Degree  of 
polynomial 

Rms  error 

Mean  absolute 

Extreme  error  (%) 

Parameter 

(XlO-6) 

error  (%) 

Minimum 

Maximum 

4 

5.2 

1.67 

-4.3 

+9.3 

7 

4.2 

1.02 

-2.4 

+4.2 

Density 

-2.4 

+4.2 

8 

3.9 

0.96 

9 

3.6 

1.02 

-2.2 

t'ij  ,9 

3 

7.6 

4.43 

-9.1 

+18.1 

5 

1.4 

0.76 

-2.8 

■ri.5 

Viscosity 

-2.9 

+1.6 

7 

1.4 

0.76 

9 

1.4 

0.74 

-3.0 

rl.6 

Because  the  errors  associated  with  a  seventh-degree  polynomial  for  density  and 
a  fifth-degree  polynomial  for  viscosity  are  of  the  same  order  of  magnitude  as  depar¬ 
tures  from  standard  atmosphere  values  caused  by  daily  changes  in  atmospheric 
structure,  and  because  little  reduction  in  error  can  be  realized  without  a  large  increase 
in  degree,  the  seventh-  and  fifth-degree  equations  were  selected.  The  percent  errors 
at  each  data  point  from  the  surface  to  89.5  km  are  shown  in  Fig.  2-1  for  density  and 
viscosity.  In  both  cases,  the  largest  error  occurred  at  the  tropopause  level  (near 
10  km)  in  the  standard  atmosphere.  The  error  profile  for  mean  free  path  is  a  mirror 
image  of  that  for  density.  Greater  accuracy  can  be  achieved  with  lower  order  poly¬ 
nomials  if  the  total  layer  is  subdivided  into  smaller  layers.  However,  this  is  not 
acceptable  in  the  present  problem  because  serious  errors  may  be  introduced  in  the 
trajectory  computations  due  to  discontinuities  in  the  time  derivatives  of  particle  iull 
speed. 
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standard  atmosphere  densities  and  viscosities^ 


At  the  user’s  option,  alternative  tabulations  of  density  and  viscosity  not  too  dif¬ 
ferent  from  the  U.S.  standard  atmosphere  values  can  be  used  in  the  model.  If  this 
option  is  used,  new  coefficients  are  computed  for  Eqs.  (2-1),  (2-2),  and  (2-3)  by  a 
least  squares  fit  to  the  desired  values. 


2.2.2  Drag-force  Equations 

The  resistance  or  drag  force  acting  on  a  particle  moving  steadily  through  the 
atmosphere  depends  on  the  Reynolds  number,  R  =  a  particle  shape  factor 

G  Tj 

B  =  the  Knudsen  number  K  =  A./d;  and;  the  Mach  number  M  =  W/v  .  In  these 
2V  n  s 

expressions,  W  is  the  particle  speed,  p  is  air  density,  d  is  particle  diameter,  77  is  air 
viscosity,  S  is  the  particle  cross-sectional  area,  V  its  volume,  X  the  mean  free  path. 


and  V  the  speed  of  sound.  For  solid  spheres,  the  shape  factor  B  has  a  minimum 
s 

numerical  value  of  3.  Following  Opik  [14]  and  Davies  [3],  the  drag  force  F  can  be 
formulated  as  follows  for  various  ranges  of  the  nondimensional  parameters: 


^2  = 


^4  = 


^3  = 


^  p(W^  +  1.5u^ 

(M  >  1; 

>  1) 

(2-4) 

(M  >  1; 

(2-5) 

yn  2  2 

n  h  R 

8p  e 

(M  >  1; 

1;  F^ 

(2-6) 

2 -5^.2 
^  d  pp  W 

(M  S  1;  K 

>Fp 

(2-7) 

^  d^  pW^ 

(M  £  1;  K 

>  1;  F^ 

(2-8) 

f /r 

8p  e 

(M  S  1;  K 
'  n 

s  1-  F 
’  3 

>F2) 

(2-9) 

VTT  2  2 

^d  pW 

(M  1; 

<  1-  F 
■  3 

"  'V- 

(2-10) 

In  the  formula  for  F  ,  the  quantity  u  is  the  root -mean-square  thermal  velocity  com- 

^  1/2  . 

ponent  along  one  coordinate  (u  =  [kT/m]  )  where  k  is  the  Boltzman  constant,  T  is 


temperature,  and  m  is  the  mass  of  one  mole  of  air.  The  parameter  y  in  the  formula 
for  Fg  is  the  drag  coefficient.  Values  of  y  recommended  by  Opik  [14]  are  given  in 
Table  2-2. 
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TABLE  2-2 

DRAG  COEFFICIENT  FOR  SHIELDED  HIGH-SPEED  FLOW  114] 


0.75/K 

n 

Mach  number 
(M) 

Drag  coefficient 
(X) 

0 

M  »  1 

2.00 

1 

M  »  1 

1.50 

2 

M  »  1 

1.25 

4 

M  »  1 

1.12 

OO 

M  »  1 

1.00 

oo 

M  1 

0.50 

Values  of  x  required  in  the  formula  for  F  are  implicit  in  Davies  interpolation  equa- 

o 

tions  for  R  [31. 
e 

2 

vRe  -4  2  2  -6  2,3 

~  ^  ^^^e  ^  ^  ^^^e  ^ 

-9  2  4 

-  6.9105  X  10  (xRg  )  (Rg  <  4) 

2  -2  2  2 

log  Rg  =  -  1.29536  +  0.986  (log  x\  )"  4.6677  X  10  (log  x\  ) 

-3  2  3  4  (2-12) 
+  1.1235  X  10  (log  xRg  )  (3  <  Rg  <  10  ) 

The  drag  coefficients  given  by  Eqs.  (2-11)  and  (2-12)  are  used  when  the  viscosity 
of  the  air  is  important,  while  those  in  Table  1  are  for  high-speed,  low-density  flow. 
When  the  Knudsen  number  is  near  unity,  Eqs.  (2-6)  and  (2-9)  can  be  used  only  if  the 
Reynolds  number  is  small,  because  the  slip  flow  corrections  for  these  equations  are 
not  known  for  large  Reynolds  numbers.  Experimental  measurements  of  the  drag 
coefficients  of  spheres  moving  at  high  speed  in  a  low-density  medium  are  being 
sought  in  an  attempt  to  check  the  theoretical  values  in  Table  2-2. 

2.2.3  Terminal  Velocities 

A  falling  particle  reaches  its  terminal  velocity  in  the  atmosphere  when  the  drag 
force  equals  the  attractive  force  of  gravity.  For  a  sphere  of  large  density  in  air, 

F  =  l  d%pg  (2-13) 
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where  p  is  the  particle  density  and  g  is  the  acceleration  of  gravity.  The  terminal 
P 

velocity  W  -  may  be  derived  from  Eq.  (2-13)  using  Eqs.  (2-4)  to  (2-12)  for  F. 
Computations  of  have  been  made  for  selected  particle  sizes  in  atmospheric  layers 
for  which  the  Mach  number  and  Knudsen  number  are  near  unity,  i.e.,  in  regions  of  the 
atmosphere  for  which  formulas  F^,  F^,  F^,  and  F^  overlap.  Terminal  velocities 
according  to  Eqs.  (2-4),  (2-7),  and  (2-11)  for  a  500  micron  diameter  sphere  (density  5) 
falling  between  60  and  90  Ion  under  standard  atmosphere  conditions  are  shov/n  in  Fig. 
2-2.  For  this  particle  size  and  altitude  range,  the  Reynolds  number  is  less  than  one, 
the  mean  free  path  exceeds  the  particle  size  above  about  59  km,  and  the  Mach  number 
exceeds  one  above  about  72  Ian.  Equation  (2-11)  yields  the  largest  drag  force  up  to 
about  78  Ian,  because  Eq.  (2-7)  is  invalid  below  59  km.  Above  78  km,  Eq.  (2-4)  for 
unshielded  supersonic  flow  is  more  appropriate  than  Eq.  (2-11).  It  is  clear  from  the 
divergence  of  the  curves  above  78  km  in  Fig.  2-2  that  there  exists  a  region  for  large 
particles  and  high  altitudes  in  which  the  terminal  velocity  equations  used  in  fallout 
Model  A  [Eqs.  (2-11)  and  (2-12)]  are  invalid.  Two  other  regions  in  which  these  equa¬ 
tions  are  invalid  have  been  identified.  The  first  of  these  first  appears  at  altitudes  of 
about  40  km  as  particle  size  is  increased  and  is  associated  with  the  conditions 

R  >  1,  K  >1.  The  second  of  these  regions  occurs  at  lower  altitudes  for 
e  n 

R  >  1,  K  <1,  and  Mach  numbers  of  the  order  of  1  or  larger.  The  validity  of  Eqs. 
e  n 

(2-5),  (2-6),  (2-8),  (2-9),  and  (2-10)  for  the  computation  of  terminal  velocity  in  these 
regions  is  under  study. 

The  behavior  of  terminal  velocities  computed  from  Eqs.  (2-7)  and  (2-11)  for  a 

smaller  particle  of  11  microns  diameter  (density  5)  is  shovm  in  Fig.  2-3.  In  this  case 

R  «  1  and  K  =  1  at  about  33  km.  Equation  (2-11)  is  selected  here  because  it 
e  n 

yields  the  greatest  drag  force  at  all  altitudes  between  33  and  90  km.  The  curves  above 
70  km  (not  shown  in  Fig.  2-3)  are  nearly  parallel. 

2.2.4  Particle  Acceleration 

Although  the  transient  acceleration  phase  experienced  by  a  particle  with  an 
initial  velocity  different  from  its  terminal  velocity  is  negligible,  in  the  lowest  layers 
of  the  atmosphere, for  most  particles  of  interest  in  the  present  problem,  it  is  sig¬ 
nificant  at  high  altitudes.  The  drag  forces  for  two  limiting  cases  of  high  altitude  flow 
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Fig.  2-2.  Terminal  velocities  based  on  dragforce  formulation  for  shielded  viscous  flow  [•,  Eq  (2-11)], 
unshielded  subsonic  flow  [a  ,  Eq  (2-7)],  and  unshielded  supersonic  flow  (o  ,  Eq  (2-4)]  for  a  500-micron  diameter 
sphere  (density  5). 


ui>i  ‘apnTOiv 


subsonic  flow  [a,  Eq  (2-7)]  for  a  11- 


are  given  by  Eq.  (2-4)  (large  W)  and  Eq.  (2-7)  (small  W).  In  the  latter  case,  the 
resistance  force  is  proportional  to  W  and  the  particle  acceleration  is  given  by 


dW 

dt 


g 


k^W 

m 


(2-14) 


where  m  is  the  particle  mass  and  is  a  constant.  Assuming  constant  atmospheric 
parameters,  this  may  be  integrated  to  give 

t/W, 


W  =  W.J,  {1  -  |1  -  w^/w.^1  exp 


(g  T) 


(2- 15) 


where  W  =  at  t  -  0.  If  the  particle  falls  from  rest,  the  characteristic  time  or 
time  constant  of  the  acceleration  phase  is  W,j,/g.  However,  with  a  non-zero  initial 
velocity,  W  the  characteristic  time  for  the  acceleration  phase  depends  also  on  the 
ratio  Solving  for  t  in  Eq.  (2-15),  we  find 


t  =  In 


(1  -  W^/W.j.)] 

(1  -  W/W.J,) 


(2-16) 


The  times  t  required  to  reach  fall  speeds  W  within  10  percent  of  the  particle  terminal 
velocity  W.^  are  shown  in  Fig.  2-4  for  two  extreme  initial  velocities  (W„  =  0  and 


0 


Wq  =  18,000  mph). 


According  to  Eq.  (2-4)  for  large  values  of  W  (i.e.,  W  >  u),  the  resistance  force 
2 


is  proportional  to  W  and  the  particle  acceleration  is  given  by 

.2 


dW 

dt 


g  - 


1^ 

m 


(2-17) 


Again  assuming  constant  atmospheric  parameters,  Eq.  (2-17)  yields,  after 
integration: 


W  =  W.J,  (exp  +  1)  ^ 


1  \V  \ 

T 

0 

-  w 

\  T 

0/ 

exp  (2gt/W.^) 


(2-18) 


where  W  =  at  t  =  0  (W^  W.^,).  Solving  for  t  in  Eq.  (2-18)  gives 


W„ 


2g 


In 


'(1  +  Wp/W,j,)  (1  -  W^/W.j,)^ 

(1  -  W/W.^)  (1  +  W^/W.^)'^ 


(2-19) 


The  times  t  required  to  reach  fall  speeds  W  within  10  percent  of  the  particle 
terminal  velocity  are  shown  in  Fig.  2-5  for  the  same  initial  velocities  as  those 


11 


Terminal  velocity  (W^),  cm  sec 


Fie  2-4  Time  required  to  attain  fall  speed  within  10  percent  of  terminal  velocity 
for  two  extreme  initial  velocities  with  drag  force  proportional  to  particle  speed. 


used  in  Fig.  2-4. 

3  -1 

Because  terminal  velocities  in  excess  of  10  cm  sec  are  found  even  for 
particles  of  the  order  of  10  microns  diameter  at  the  highest  altitudes  of  interest,  it  is 
apparent  from  Figs.  2-4  and  2-5  that  rather  severe  limitations  must  be  placed  either 
on  particle  size  or  altitude  unless  the  acceleration  phase  is  taken  into  account. 

Although  the  characteristic  times  in  Figs.  2-4  and  2-5  indicate  that  particle 
accelerations  are  important,  these  accelerations  cannot  be  used  to  estimate  the  effect 
of  acceleration  on  particle  residence  time  because,  due  to  changes  in  air  density,  vis¬ 
cosity,  and  mean  free  path,  the  particle  terminal  velocities  decrease  rapidly  with  a 
decrease  in  altitude.  Numerical  integrations  of  the  fall-speed  equations,  which  take  into 
account  both  the  accelerations  and  the  changes  in  drag  forces  with  altitude,  have  been 
undertaken  in  an  effort  to  identify  the  upper  limits  of  particle  size  and  altitude  for  which 
acceleratiohs  may  be  neglected  in  terms  of  effect  on  total  residence  time  or  net  hori¬ 
zontal  displacement  in  deep  layers.  In  the  reentry  fallout  problem,  the  particle 
accelerations  may  be  of  importance  only  if  they  lead  to  appreciable  changes  in  the 
time  of  arrival  or  position  of  particles  at  levels  of  10  km  or  lower. 

2.3  Mass-element  Trajectories 

2.3.1  Introduction 

The  positions  of  the  mass  elements  or  subdivisions  of  each  particle  source- 
cloud  element  must  be  known  at  specified  levels  in  the  atmosphere  to  compute  air 
concentrations  and  deposit  densities  at  the  earths’  surface.  The  equations  for  the 
position  of  the  center  of  a  mass-element  at  the  (r  +  1)  time  step  are: 


r+1 


0  + 


f^r 

T 

r 


0(0,  z,  t)  dt 


(2-20) 


^  1  ^  A.  (0,  A,  z,  t)  dt 

r+1  r  t 

r 


(2-21) 


z^^^  =  z^  +  J^r  ^  {  W(z)  +  w(9.  A,  z,  t)}dt  (2-22) 

r+  r  tj, 

where  9  (colatitude),  A  (longitude),  and  z  (altitude)  are  position  coordinates  in  the 
atmosphere,  t  is  time,  W  (z)  is  the  terminal  velocity  of  particles  in  a  mass-element. 


14 


0  and  X  are  the  meridional  and  zonal  components,  respectively,  of  the  an^idar  velocity 
of  the  air,  and  w  is  the  vei’tical  air  velocity.  It  is  assumed  that  the  horizontal  particle 
velocity  is  equal  at  all  times  to  the  local  horizontal  air  velocity  and  that  the  vertical 
particle  velocity  is  equal  to  the  sum  of  the  particle  terminal  velocity  and  the  local 
vertical  air  velocity.  The  latter  assumption  was  discussed  in  Section  2.2,4. 

Because  of  the  complexity  of  the  analytic  forms  for  0,  X,  W,  and  w,  it  is  neces¬ 
sary  to  evaluate  Eqs.  (2-21),  (2-22),  and  (2-23)  by  numerical  methods.  Numerous 
methods  of  integration  are  available  and  it  becomes  a  matter  of  selecting  the  optimum 
one  for  this  problem. 

The  principal  selection  criteria  are; 

(a)  The  position  of  the  end-point  of  the  trajectory  must  be  insensitive  to 
such  details  of  the  method  of  computation  as  the  choice  of  time-step  length, 
the  choice  of  criteria  for  the  control  of  individual  time  step  errors,  etc. 

(b)  It  should  be  possible  to  estimate  reliably  and  with  case  the  initial 
time-step  length  and  subsequent  changes  in  time-step  length  as  required  to 
achieve  (a). 

(c)  The  eomputation  time  should  be  minimized. 

(d)  The  programming  and  verifying  procedures  should  be  as  straight¬ 
forward  and  simple  as  possible, consistent  with  (a),  (b),  and  (c). 

The  accuracy  requirements  described  by  (a)  above  can  be  determined  quanti¬ 
tatively  by  comparing  the  accumulated  position  errors  with  estimates  of  uncertainty- 
due  to  natural  diffusion  processes.  If  the  particle  position  errors  due  to  computation 
are  large,  an  artificial  particle-cloud  growth  rate  will  result  and  the  estimates  of 
concentrations  and  deposit  densities  will  be  meaningless. 

Two  specific  methods  of  integration  representing  the  one-step  Runge-Kutta 
method  and  the  multi-step  predictor-corrector  class  of  methods  were  selected  for 
testing.  Both  methods  were  programmed,  checked,  and  subjected  to  preliminary 
performance  tests  on  simplified  equations  possessing  some  of  the  properties  of  the 
equations  that  will  later  be  used  to  simulate  the  atmospheric  wind  field.  This  work  is 
described  in  more  detail  in  the  following  sections. 
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2.3.2  Runge— Kutta  Integration  Method 

Equation  (2-22)  for  particle  height  will  be  used  to  illustrate  the  Runge— Kutta 
procedure.  If  h  is  the  time  step, 


z(t^  +  h)  =  z{t^)  +  W[t,z(t)]  dt  (2-23) 

omitting  the  term  w  (0,  A,  z,  t)  for  simplicity.  The  Runge-Kutta  formulas  represent 
an  approximation  (to  a  specified  order  of  accuracy)  to  a  formal  Taylor’s  series  expan¬ 
sion  for  z  (t^  +  h).  The  approximation  requires  the  evaluation  of  first  derivatives 
only.  The  following  quantities  are  needed  for  a  fourth  order  approximation,  i.e.,  exact 

4 

agreement  with  a  Taylor  series  expansion  through  terms  of  order  h  : 


a^  =  hW(t^.  z^) 

Ea  =  hW(tQ  +  b^h,  Zq  +  c^a^  +  d^a^) 

\  =  hW(tQ  +  bgh,  +  Cga^  -  d^a^  e^Ea) 

Z(to  +  h)  -  =  f^a^  .  f^a^  -  f3aa  .  f^a^. 


(2-24) 


(2-25) 


The  form  of  Eqs.  (2-24)  and  (2-25)  represents  a  generalization  of  the  form  that  results 
from  the  approximation  of  a  Taylor’s  series  solution  by  Simpson’s  rule.  The  symbol 
Z  indicates  that  this  quantity  is  an  approximation  to  z  (t^  +  h). 

All  but  one  of  the  parameters  b,  c,  d,  e,  and  f  are  needed  to  meet  the  require  ¬ 
ment 

Z  (t^  +  h)  =  Z  (t^  +  h)  +  0  (h"^) . 

Gill  [8],  in  optimizing  the  method  for  use  on  a  computer,  chose  the  arbitrary  condition 
to  be  e^  =  1  +  /1/2.  The  other  parameters  are  then 

b^  =  1/2  ^ 

b„  =  1/2  ^^9  =  - 

=  1/2  fg  =  1/3  (1  -  Vl72') 
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Cg  =  -  1/2  +  vT:^2  =  1/3  (1  +  vT/2) 

C3  =  0  =  1/6 

The  form  of  the  Runge— Kutta— Gill  equations  used  in  the  computer  program  involves 
auxiliary  quantities,  g,  as  recommended  by  Gill  [8]  for  greater  efficiency  and  round¬ 
off  error  control.  The  solution  is 

=  Z3  +  1/6  (a^  -  2g3)  =  Z  (t^  +  h)  (2-26) 

where 


a 

1 


a 


3 


a 


4 


=  hW(t^,  z^) 

=  hW(z^  +  h/2,  z^) 
=  hW(z^  +  h/2,  Zg) 
=  hW(z^  +  h,  Zg) 


z  +  1/2 (a  -  2g  ) 
0  1  0^ 

z^  +  (1  -  im)  (a^ 

z^  +  (1  +  fl72)  (a^ 


gg) 


(2-27) 


=  gp  +  3  |l/2  (a^  -  2g^)]  -  1/2 a^ 

gg  =  g^  +  3  [(1  -  vT/2)  (a^  -  g^)]  -  (1  -  /172)  a^ 

gg  ="  gg  +  3  [(1  +  %'l72)  (Ug  -  g^)]  -  (1  +  Vl72)ag 

"  *^3  ^  *^4  "  ^®3^^  “  ^4 

with  g^  =  0  in  the  first  time  step  and  g^  =  g^  in  subsequent  time  steps. 

The  Runge-Kutta  method  described  above  was  coded  for  a  system  of  three 
equations  (9,  X,  and  z)  following  the  procedure  described  by  Ralston  and  Wilf  [15]. 


The  coding  for  the  9  and  X  equations  was  checked  using  a  sample  problem  for  which 


the  solutions  were  readily  available  [15,  p.  119].  Verification  of  the  coding  for  the 
solution  of  the  z  equation  is  described  in  a  later  section. 


2.3.3  Predictor— Corrector  Integration  Method 

Predictor-corrector  techniques  for  the  solution  of  an  ordinary  differential  equa- 
dz 

tion  such  as  ~  =  W(z,t)  involve  the  prediction  of  a  value  of  z  at  t  =  t^  +  h,  the 

dz 

computation  of  an  estimate  of  ~  at  t  =  t^  +  h  using  the  predicted  z,  and  the  com¬ 
putation  of  a  corrected  approximation  to  z  at  t  =  t^  +  h.  Hamming’s  method  [10] 
which  is  both  accurate  and  stable  for  many  integration  problems  was  coded  for  testing 
on  the  trajectory  problem.  In  this  method,  the  corrector  is  derived  by  evaluating  the 
coefficients  of  a  linear  combination  of  z  (t^),  z(t^  -  h),  zjt^  -  2h),  z(t^  +  h),  zjt^), 
and  z(tQ  -  h)  from  a  Taylor’s  series  expansion  of  z(tQ  +  h)  under  the  condition  that 
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all  powers  of  h  through  h'^  vanish.  The  accuracy  of  the  predictor  is  greatly  increased 
by  the  use  of  a  modifier  derived  from  an  estimate  of  the  truncation  error  at  each  step. 
Using  the  notation  p  for  predictor,  m  for  modifier,  and  c  for  corrector,  and  numerical 
subscripts  to  identify  the  time  of  evaluation  of  each  term  relative  to  t^,  the  integration 


equations  are: 


Pi  =  "-3 


4h 


(2z„  -  z  -I-  2z  „) 
3  '  0  -1  -2' 


112 


m. 


Pi  121  ^Pq 


m^  =  W(t^,mp 


h  ■  i  (‘'"o  "  "-2  *  -  "-!» 


"i  “  h  *  lii  ">1  ■ 


(2-28) 


(2-29) 


The  above  method  has  an  important  advantage  over  the  Runge-Kutta  method 
because  it  provides  an  estimate  of  the  error  at  each  step  and  these  estimates  can  be 
used  to  change  the  step  length.  However,  the  method  is  not  self-starting.  In  the  first 
version  that  was  programmed,  the  Runge-Kutta  method  was  used  exclusively  for  both 
starting  and  changing  step  length.  However,  due  to  the  need  for  frequent  changes  in 
step  length  in  test  trajectory  problems,  this  approach  did  not  provide  a  suitable  basis 
for  comparing  the  performance  of  the  two  integration  techniques.  Consequently,  a 
second  version  using  stored  values  of  z  for  doubling  the  time  step  and  interpolated 
values  of  z  for  halving  the  time  step  was  prepared.  The  interpolation  formulas  quoted 
by  Ralston  and  Wilf  [15,  p.  103]  were  used  because  they  possess  the  same  order  of 
accuracy  as  the  integration  method. 


-  1/2  =  ~  (80  +  135  z_^  +  40  z_2  +  Z_g) 


—  (-15  z^  +  90  z_^  +  15  z_2) 
1 


(2-30) 


-3/2  256 


(12  z^  +  135  z_^  +  108  z_2 


+  —  (-3  z 
256  '■  n 


54  z  +  27  z  ) 
n-1  n-2 


The  coding  of  the  predictor-corrector  integration 
sample  problem  given  by  Ralston  and  Wilf  [15,  p.  108]. 


method  was  checked  using  the 
Preliminary  comparisons  of 
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the  performance  of  the  Runge-Kutta  and  predictor-corrector  methods  arc  described 
in  the  following  section. 

2.3.4  Integrations  of  Some  Simplified  Particle  Trajectories 

Tests  of  the  two  integration  methods  were  performed  first  on  Eq.  (2-23)  using 

Davies  equations  [3]  for  particle  terminal  velocity,  W[t,z(t)],  and  Eqs.  (2-1),  (2-2), 

and  (2-3)  for  the  vertical  distribution  of  density,  viscosity,  and  mean  free  path. 

Typical  results  are  shown  in  Table  2-3,  which  lists  computed  heights  as  a  function  of 

time  for  a  97-micron  diameter  sphere  (density  5)  falling  from  29  km.  No  changes  in 

step  length  were  permitted  during  integration  by  the  predictor-corrector  method  in 

this  test.  The  first  few  values  are  identical  for  both  methods,  because  the  Runge- 

Kutta  method  was  used  for  starting.  It  is  clear  from  Table  2-3  that  both  methods 

4 

yielded  accurate  results  even  for  the  largest  time  step.  At  t  =  2  x  10  sec,  the 

greatest  variation  in  height  estimates  was  about  0.05  km  or  less  than  0.2  percent  of 

4 

the  total  distance  of  fall.  The  fall  time  of  2.10  x  10  sec  from  29  km  to  0.3  km  com- 

4 

pares  reasonably  well  with  the  total  time  of  2.19  x  10  sec  obtained  with  Model  A  for 
the  fall  of  a  similar  sphere  from  30  Ion  to  the  earth’s  surface. 

For  most  particle  sizes  of  interest  the  terminal  velocity  increases  monotonically 
with  altitude.  The  zonal  wind  speed  k,  on  the  other  hand,  may  be  characterized  by 
several  minima  and  maxima  representing  easterly  and  westerly  wind  layers.  The 
following  equation  was  used  to  represent  X  for  initial  test  purposes: 

X  =  A  sin  kz  (2-31) 

6  “1  1 

with  the  maximum  wind  speed  A  =  7.87  x  10  rad  sec  (50  m  sec  )  and 
k  =  2n/15.  The  wind  profile  given  by  Eq.  (2-31)  is  shown  in  Fig.  2-6  (for  latitude  0“). 

Using  Davies’  [3]  equations  for  terminal  velocity  and  Eq.  (2-31)  for  zonal  wind 
speed,  the  height  z  and  longitude  X  were  computed  for  a  97-micron  diameter  sphere 
(density  5)  as  functions  of  time  by  the  Runge-Kutta  and  predictor-corrector  methods. 

3 

The  results  for  time  step  h  =  10  sec  are  listed  in  Table  2-4.  The  criteria  used  for 
changing  time-step  length  in  the  predictor-corrector  method  were  as  follows: 

2  /  1p^  -  cM  m  (i  =  1,  2,  3)  (2-32) 

i  r  r  1 
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TABLE  2-3 

PARTICLE  ALTITUDE  [z,  (km)]  AS  A  FUNCTION  OF  TIME  (t)  AS  COMPUTED 
BY  NUMERICAL  INTEGRATION  USING  THE  RUNGE-KUTTA  (R-K)  AND 
PREDICTOR-CORRECTOR  (P-C)  FOURTH-ORDER  METHODS  WITH  VARIOUS 


TIME  STEPS  (h).  (97-micron  diameter  sphere,  density  5) 


t 

h  = 

100 

h  = 

500 

h  = 

1000 

h  = 

5000 

(X  10“^  sec) 

R-K 

P-C 

R-K 

1 

n 

R-K 

P-C 

R-K 

P-C 

0 

29.000 

29.000 

29.000 

29.000 

29.000 

29.000 

29.000 

29.000 

1 

27.211 

27.211 

27.211 

27.211 

27.211 

27.211 

- 

- 

2 

25.452 

25.452 

25.452 

25.452 

25.452 

25.452 

- 

- 

3 

23.729 

23.729 

23.729 

23.729 

23.729 

23.729 

- 

4 

22.046 

22.046 

22.046 

22.046 

22.046 

22.046 

- 

- 

5 

20.403 

20.404 

20.403 

20.407 

20.408 

20.404 

20.392 

20.392 

6 

18.797 

18.797 

18.796 

18.800 

18.801 

18.793 

- 

- 

7 

17.239 

17.240 

17.239 

17.242 

17.244 

17.238 

- 

— 

8 

15.730 

15.731 

15.730 

15.733 

15.734 

15.728 

- 

— 

9 

14.267 

14.267 

14.266 

14.269 

14.271 

14.262 

- 

— 

10 

12.843 

12.844 

12.842 

12.845 

12.847 

12.841 

12.833 

12.833 

11 

11.447 

11.448 

11.445 

11.454 

11.456 

11.448 

- 

- 

12 

10.100 

10.101 

10.099 

10.106 

10.108 

10.093 

- 

13 

8.8182 

8.8194 

8.8170 

8.8239 

8.8264 

8.8185 

- 

— 

14 

7.6001 

7.6013 

7.5990 

7.6055 

7.6079 

7.5982 

- 

- 

15 

6.4410 

6.4420 

6.4398 

6.4461 

6.4482 

6.4408 

6.4401 

6.4401 

16 

5.3275 

5.3286 

5.3264 

5.3324 

5.3345 

5.3264 

- 

— 

17 

4.2552 

4.2562 

4.2541 

4.2599 

4.2619 

4.2543 

- 

18 

3.2201 

3.2211 

3.2191 

3.2246 

3.2266 

3.2191 

- 

— 

19 

2.2190 

2.2199 

2.2180 

2.2234 

2.2253 

2.2181 

- 

— 

20 

1.2491 

1.2500 

1.2481 

1.2533 

1.2552 

1.2482 

1.2483 

1.2860 

21 

0.3079 

0.3088 

0.3070 

0.3120 

0.3138 

0.3070 

— 

— 
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?  f 

1 

1  i  i 

Ip  -  c 
r  r 

Um2 

(i  = 

1,  2,  3) 

(2-lVS) 

? 

1  i  i  1 

Ip  -  c  i 

*•  v*  v» 

>m„ 

(i  = 

1,  2,  3) 

(2-34) 

where  |  -  cM  represents  the  absolute  value  of  the  difference  between  the  predictor 

and  the  corrector  at  the  r^^  time  step  along  coordinate  i  and  where  f  and  m  are  con¬ 
stants  assigned  the  following  numerical  values; 


/  (z)  =  1.0 
/  (X)  =  0.08 
0 


m. 


m. 


6x10 


6  X  10 


-4 


3 

i  (9) 


m  =  3  X  10 

O 


The  constants  £  above  are  weighting  factors  that  permit  proper  apportionment  of  the 
individual  step  errors  among  the  three  coordinates.  The  constants  m^,  m^,  and  m^ 
represent  the  criteria  for  doubling  the  time  step,  halving  the  time  step,  and  stopping 
computations,  respectively. 

3 

Time  step  length  changes  occurred  at  t  =  4  x  10  sec  (doubled)  and  again 

4 

at  t  =  10  sec  (halved)  with  the  predictor— corrector  method  in  Table  2-4.  As  expected, 
the  oscillatory  nature  of  the  zonal  wind  profile  resulted  in  much  larger  differences 
between  simultaneous  values  of  X  computed  by  the  two  methods  than  those  found  for  z 
(see  Table  2-3).  The  zonal  positions  of  the  particle  given  by  the  Runge-Kutta  integra¬ 
tions  are  shown  in  Fig.  2-6  (right  side). 

Because  the  integrations  listed  in  Table  2-4  are  based  on  Davies’  equations  for 
terminal  velocity,  the  correct  solutions  for  z  and  X  are  not  known  exactly.  Thus,  the 
following  simplified  equation  was  used  for  terminal  velocity; 


i/t 


(2-35) 


where  n  is  a  constant.  Equations  (2-31)  and  (2-35)  can  be  integrated  to  give  the  follow¬ 
ing  explicit  expressions  for  z  and  X; 


z  =  Q  In  t  +  C. 


X  = 


At 


[1  +  (k  Q)  ] 
AkOt 
[1  +  (kQ)^] 


cos  k  {sin  (k  Qfn  t)  -  (kQ)  cos  (kQfn  t)} 


(2-36) 


(2-37) 


sin  k  {sin  (kQjfnt) 


+  — :  cos  (k Qfn t)) 

(kQ)  ^ 


C. 
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Fig.  2-6.  Idealized  zonal  wind  profile  (left)  and  particle  trajectory  (right)  computed  by  Runge — Kutta 
numerical  integration  (time  step  =  10"^  sec). 


where 


z  -  Q  M 
0  0 


(2-38) 


and 


C 


1 


s  = 


- -  {sinkz^  -  (kO)  cos  kz^} 

[1  +  (kQ)  ] 


(2-39) 


The  various  constants  in  Eqs.  (2-36)  and  (2-37)  can  be  interpreted  as  follows. 

A  particle  is  released  at  time  t^  from  altitude  z^  with  an  initial  terminal  velocity  0,. 
It  falls  at  decreasing  speed  through  alternating  layers  of  easterly  and  westerly  winds. 
The  maximum  zonal  wind  speed  in  each  layer  is  given  by  A  and  the  total  number  of 
wind  layers  between  the  earth’s  surface  and  z^  is  kz^/n.  For  the  choice  z^  -  30  km, 
t^  =  1.0  hr,  and  Q=  -30  km  hr”^,  the  terminal  velocity  varies  with  altitude  according 
to 


W„ 


30  exp 


30 


(2-40) 


Equation  (2-40)  is  illustrated  in  Fig.  2-7.  The  remaining  constants  were  assigned  the 
values  A  =  4.68  x  lo'^  rad  hr’^  and  k  =  2n/15.  The  error  criteria  [Eqs.  (2-32), 
(2-33),  and  (2-34)]  were  chosen  to  be: 


f  (z)  =  1.0 

2  2 
I  (\)  =  6.4  X  10 

(0)  =  0 


m^  =  0.003 
m^  =  0.03 
m  =  0.15 

O 


The  zonal  components  of  the  particle  trajectory  computed  by  the  Runge-Kutta 
and  predictor— corrector  methods  using  four  different  initial  time  steps  are  given  in 
Table  2-5.  The  correct  longitudes  (rounded  to  5  significant  figures)  from  Eq.  (2-37) 
are  listed  in  the  last  column  of  the  table.  Due  to  changes  in  step  length  during  the 
integration,  predictor— corrector  estimates  were  not  available  for  all  selected  times. 
The  number  of  time  steps  in  the  bottom  row  of  the  table  refers  to  all  computations 
between  t  =  1.000  and  t  =  2.624  hr. 

For  the  largest  time  step  (h  =  0.232  hr),  predictor— corrector  computations 
were  terminated  at  t  =  2.1  hr  due  to  excessive  errors.  In  terms  of  both  accuracy  and 
speed,  the  performance  of  the  Runge-Kutta  method  was  superior  to  the  predictor 
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Number 
of  time 
steps 


correctoi’  method  for  this  problem.  In  computations  performed  by  the  latter  method, 
large  errors  appeared  first  during  changes  in  step  length. 

For  a  final  series  of  trajectory  integration  tests,  a  third  equation  for  north- 
south  particle  motions  was  added  to  those  previously  used  for  the  east-west  and 
vertical  components  of  motion. 

e  =  B  sin  Kt  sin  MG  (2-41) 

Integration  of  Eq.  (2-41)  gives 


where 


tan 


exp  -  M 


cos  Kt 


(2-42) 


M 


MG 


in  tan  — 


B 

K 


cos  K 


(2-43) 


The  constant  B  is  the  maximum  meridional  wind  speed.  The  number  of  complete  waves 
in  the  north- south  component  of  particle  motion  between  the  initial  altitude  and  the 
earth’s  surface  is  Kz^/u. 

A  total  of  four  tests  were  performed  using  the  Runge-Kutta  method  with  variations 
in  0Q,  B,  and  the  step  length  h.  Input  parameters  for  the  tests  are  listed  in  Table  2-6. 
The  values  of  k  and  K  in  Table  2-6  specify  six  wind  layers  in  the  vertical  and  twenty 
waves  in  the  meridional  component  of  motion  of  the  particle  between  90  km  and  the 
earth’s  surface.  The  initial  co-latitude,  0^,  and  maximum  meridional  wind  speed,  B, 
in  Tests  1  and  2  were  chosen  to  yield  large  round-off  errors  in  order  to  test  the 
behavior  of  the  integration  methods  under  such  conditions.  The  results  of  Tests  1  and 
2  are  given  in  Table  2-7. 

It  is  evident  from  Table  2-7  that  the  numerical  integrations  were  less  accurate 
for  the  oscillatory  functions  for  G  and  X  than  for  the  monotonic  function  for  terminal 
velocity  The  errors  in  X  appeared  to  be  due  primarily  to  truncation,  because 
their  magnitudes  increased  as  the  time  step  length  was  increased.  The  opposite 
occurred  for  0  errors,  however,  as  they  were  due  primarily  to  round-off.  The  Runge- 
Kutta  estimates  of  G  were  well-behaved,  with  no  indication  for  instability  even  with  a 
relatively  large  time  step  length.  The  time  step  length  h  =  0.24  hr  corresponded  to 
approximately  0.24  x  (20/12)  =  0.4  times  the  wave-length  of  the  oscillations  in  the 
north— south  component  of  motion  of  the  particle. 
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TABLE  2-6 

INPUT  PARAMETERS  FOR  THREE-DIMENSIONAL 
TRAJECTORY  COMPUTATIONS 


In  Tests  3  and  4,  the  parameters  9^  and  B  were  increased  to  reduce  round-off 
errors.  The  results  of  exact  and  numerical  integrations  of  0  (Runge-Kutta  method) 
for  this  problem  are  listed  in  Table  2-8.  Here  the  0  errors  increased  with  an  increase 
in  step  length,  but  even  with  a  step  length  of  0.24  hr,  the  largest  errors  were  not  more 
than  3  percent  of  the  maximum  deviation  of  the  particle  from  the  source  latitude. 

Tests  1  to  4  will  be  repeated  with  the  predictor— corrector  method  of  integration 
in  December.  Following  this,  tests  will  be  performed  using  the  simulated  zonal  wind 
fields  discussed  in  Section  3.0  of  this  report. 
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PARTICLE  POSITIONS  (z,  0)  AS  A  EUNCTION  OF  TIME  (1)  BY  EXACT  INTEGRATION  AND  BY 

RUNGE-KUTTA  (R-K)  NUMERICAL  INTEGRATION  USING  TWO  TIME  STEPS  (h) 
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TABLE  2-8 

PARTICLE  COLATITUDE  (0)  AS  A  FUNCTION  OF  TIME  (t)  BY  EXACT 
INTEGRATION  AND  BY  RUNGE-KUTTA  (R-K)  NUMERICAL  INTEGRATION 

USING  TWO  TIME  STEPS  (h) 


t 

z 

0  - 

0^  (X  10“^  rad) 

(hr) 

(km) 

L 

exact 

h  =  0.12 

h  =  0.24 

0 

0 

0 

1.48 

75.9 

+  7.495 

+  7.507 

+  7.655 

1.96 

65.8 

+11.028 

+11.045 

+11.268 

2.44 

57.9 

+  7.984 

+  7.998 

+  8.155 

2.92 

51.4 

+  0.614 

+  0.622 

+  0.618 

3.40 

45.9 

-  5.629 

-  5.771 

3.88 

41.2 

-  6.132 

-  6.131 

-  6.284 

4.36 

-  0.512 

-  0.540 

4.84 

33.2 

+  7.055 

+  7.073 

+  7.212 

5.32 

29.8 

+10.993 

+11.017 

+11.239 

26.7 

+  8.378 

+  8.402 

+  8.565 

6.28 

23.8 

+  1.148 

+  1.164 

+  1.167 

6.76 

21.2 

-  5.353 

-  5.344 

-  5.484 

7.24 

18.7 

-  6.316 

-  6.305 

-  6.469 

7.72 

16.4 

-  1.010 

-  1.058 

14.2 

+  6.598 

+  6.625 

+  6.747 

8.68 

12.2 

+10.925 

+10.957 

+11.175 

9.16 

+  8.751 

+  8.781 

+  8.952 

9.64 

8.4 

+  1.685 

+  1.705 

+  1.722 

6.7 

-  5.034 

-  5.167 

-  6.469 

-  6.453 

-  6.621 

3.4 

-  1.527 

-  1.503 

-  1.563 

11.56 

1.9 

+  6.124 

+  6.160 

+  6.266 

12.04 

0.4 

+10.825 

+10.866 

+11.076 

Number  of 
time  steps 

- 

92 
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2.4  Mass-element  Expansion  and  Subdivision 

Computer  program  specifications  were  prepared  for  the  subdivision  of  mass- 
elements  forming  both  polydisperse  and  monodisperse  particle  clouds  as  outlined  in 
the  first  quarterly  progress  report  issued  under  this  project  [9].  Coding  of  this  sub¬ 
routine  was  delayed  pending  completion  of  tests  of  trajectory  integration  methods. 

2.5  Surface  Deposit  Densities  and  Air  Concentrations 

Preliminary  specifications  and  coding  were  completed  for  establishing  a  grid  at 
the  earth’s  surface  and  at  some  level  above  the  earth’s  surface  for  the  purpose  of 
analyzing  surface  deposit  densities  and  air  concentrations,  respectively.  A  polar 
stereographic  projection  is  used  for  latitudes  poleward  of  20  degrees  and  a  Mercator 
projection  for  tropical  latitudes.  Verification  of  these  components  of  the  model  have 
also  been  delayed  until  the  trajectory  routines  have  been  completed. 
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3 .0 _ ATMOSPHERIC  AIR  FLOW  SIMULATION  MODELS 

Zonal  Flow  Representation 

Examination  of  the  least-squares  polynomial  approximations  of  the  pole-to-pole 
zonal  wind  profiles  obtained  by  the  methods  described  in  [9]  has  shown  the  most 
desirable  fit  to  be  the  ninth-degree  approximation.  The  choice  of  this  polynomial  was 
determined  by  its  over-all  accuracy  in  the  latitudinal  regions  of  greatest  importance 
(70“  N  to  30“  S)  for  the  purposes  of  the  project.  In  general,  lower-degree  approxima¬ 
tions  were  found  to  give  better  results  in  the  polar  regions  at  the  expense  of  close 
agreement  with  the  observed  profiles  in  the  lower  latitudes. 

Figures  3-1  and  3-2  show  the  observed  500-  and  200-mb  profiles  for  two  sea¬ 
sons  as  given  by  Mintz  [12]  plotted  together  with  the  corresponding  9th-  and  5th-  degree 
approximations.  It  is  seen  that  an  undesirable  fluctuation  occurs  in  the  northern  polar 
region  for  the  9th-degree  curves,  but  that  the  9th-degree  polynomial  gives  an  accept¬ 
able  approximation  for  the  remaining  latitudes  while  the  5th-degree  polynomial  docs  not. 

Figures  3-3  and  3-4  show  the  approximating  profiles  for  the  30-  and  60-km 
levels.  Observed  profiles  were  not  available  for  these  upper  levels,  but  data  were 
obtained  [10]  for  the  northern  hemisphere  (winter  and  summer)  at  various  latitudes  and 
interpolated  to  equally  spaced  intervals.  Because  no  southern  hemisphere  data  of  this 
tyiDe  were  found,  data  for  the  northern  hemisphere  winter  were  taken  for  the  southern 
hemisphere  winter,  etc.  to  obtain  complete  pole-to-pole  profiles. 

For  each  of  the  specified  height  levels  z  =  z^,  (i=l,  2,  ...,  f,  where  f  is  the  total 
number  of  levels  at  which  data  are  read  off),  the  ninth -degree  approximation  to  the 
zonal  wind  has  the  form 

9 

Vg  (Z^,  0)  =  aj^(Z^)  Pj^(0)  (i“l>  •••> 

where  0  denotes  colatitude.  For  a  fixed  value  of  k,  each  set  of  coefficients  a.^ 

a  (z  ),  ...,  a  (z,  )  can  be  fitted  by  the  same  least-squares  process  as  was  used  in  the 
Ic  2  k  i 

horizontal  approximations  to  obtain  a  smooth  vertical  fitting  of  the  pole-to-pole  pro¬ 
files.  This  method  of  fitting  eliminates  the  possibility  of  encountering  undesirable 
effects  upon  particle  trajectories  brought  about  by  discontinuities  in  the  vertical  shear. 

In  order  to  take  advantage  of  the  greater  resolution  of  data  at  the  lower  altitudes. 
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V,  m  sec 

Fig.  3-1  a,  b.  The  observed  500-  and  200-mb  pole-to-pole  zonal  wind  profiles  for  the  northern-hemisphere  summer 
and  southern-hemisphere  winter,  plotted  together  with  the  5th-  and  9th-degree  approximations. 


mb  profile 
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Fig.  3-2  a,  b.  The  observed  500-  and  200-mb  pole-to-pole  zonal  wind  profiles  for  the  northern-hemisphere 
winter  and  southern-hemisphere  summer  plotted  together  with  the  5th-  and  9th -degree  approximations. 


Fig,  3-3.  The  30-km  observed  zonal  wind  values  plotted  together  with  the  5th- 
and  9th-degree  approximations. 
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[Vofile  (Kantor  ancl  Cole  i  |  ) 


observed  zonal  wind  values  plotted  together  with  the  5th-  and  9th-degree  approximations 


a  change  of  height  variable  from  z  to  z'  can  be  made  to  place  the  levels  z^,  z^, 

z,'  at  approximately  equal  spacing  in  the  new  variable  z/.  The  transformation 

^  1/2 

z'  =  z  accomplishes  this  purpose. 

A  computer  program  has  been  completed  for  calculating  the  zonal  wind  at  any 
number  of  coordinate  points  (6,  z);  thus  vertical  as  well  as  pole-to-pole  profiles  can  be 
readily  calculated  for  comparison  with  observational  data.  Analysis  of  many  such  pro¬ 
files  will  be  necessary  for  determining  the  appropriate  degree  of  approximation  in  the 
vertical  least-squares  fit. 

3.2  Long  Atmospheric  Waves 

In  accordance  with  the  formula 


T  (G,  X)  =  2  2  cos  [m(X-6^^)]  P  (cos  S) 
'  m  n  n  '  n  n 


(3-2) 


the  stream  function  for  the  long-atmospheric  waves  has  been  computed  at  e\eij  ten 
degi’ees  of  latitude  and  longitude  for  the  1000-  and  500-mb  levels.  The  amplitudes  A^ 
and  phase  angles  6^  which  were  used  in  the  calculations  are  those  given  by  Eliasen 
and  Machenhauer  [6]  for  a  90-day  winter  mean  flow.  .As  is  custoniarj ,  the  long-vave 
numbers  correspond  to  the  values  m  =  1,  2,  3,  4,  while  in  this  case  (n-m)  takes  on 
the  values  1,  3,  5,  7.  (Recall  that  m  represents  the  number  of  waves  around  a  circle 
of  latitude,  while  n-m  is  the  number  of  nodal  points  of  (9)  for  0  <  -  <  n.) 

A  computer  subroutine  for  evaluating  the  associated  Legendre  functions  P^ 

(cos  0)  was  developed  making  use  of  the  formula 


m. 


ni 

P  (cos  0) 
n 


(2n) !  sin 
2^n!  (n-m)! 


n-m  „  (n-m)(n-m-l)  n-m-2^ 

Icos  0  -  -  -  cos  0 

‘  2(2n-l) 


(n-m)  (n-m-1)  (n-m-2)  (n-m-3)  ^^^n-m-4^ 
(2)  (4)  (2n-l)  (2n-3) 


(3-3) 

where  the  expression  in  the  brackets  ends  with  the  a  constant  term  if  n-m  is  e\en  and 
with  the  term  containing  cos  9  if  n-m  is  odd.  This  equation  is  that  given  b\  Bveilv 
[2]  and  is  basically  the  same  as  the  one  recommended  for  computational  purposes  bv 
Fougere  [7].  For  verification  of  the  results,  values  of  the  normalized  associated 
Legendre  function  P^  ,  which  is  related  to  P^^  by  the  foi  mula 
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surface,  the  stream  lines  for  the  long-wave  components  are  altered  onlj-  slightly  by 
the  addition  of  '1' ;  for  this  reason,  the  map  showing  '1'  alone  has  been  omitted. 

At  the  500-mb  level,  the  effect  of  the  zonal  wind  is  substantial;  Fig.  3-6  shows 
the  stream  lines  for  the  long  waves  alone,  while  Fig.  3-7  shows  the  stream  lines  for 
the  combined  flow.  In  the  latter  case,  the  zonal  wind  dominates  the  flow  to  such  an 
extent  that  nearly  all  the  individual  cells  are  smoothed  into  ridges  and  troughs,  while 
only  two  closed  lows  are  seen  to  retain  their  identities. 

In  Fig.  3-8,  the  stream  lines  for  the  combined  flow  at  the  1000-  and  500-mb 
levels  are  shown  for  the  low’-latitude  belt  from  the  equation  to  30“  N. 

The  vertical  fitting  of  the  long-wave  components  is  to  be  done  in  a  manner 
analogous  to  that  of  the  zonal  wind  approximations,  i.e.,  and  are  to  be  taken 
as  functions  of  the  vertical  coordinate  z.  Due  to  a  lack  of  available  data  in  regard  to 
the  long-wave  components  at  the  higher  levels,  the  vertical  fitting  must  be  of  a  less 
precise  nature  than  that  described  in  the  preceding  section;  nevertheless,  some 
information  concerning  the  behavior  of  the  amplitudes  and  phase  angles  with  height  is 
available  (Muench  [13]  and  Eliasen  [5])  for  guidance  in  choosing  an  appropriate  fit. 
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Fig.  3-5,  The  stream-lines  of  the  northern-hemisphere  zonal  wind  with 
long-wave  wind  components  for  the  winter  season  at  the  1000-mb  level. 
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Fig.  3-6.  The  stream-lines  of  the  northern -hemisphere  long-wave  wind  com¬ 
ponents  for  the  winter  season  at  the  500-mb  level. 
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Fig.  3-7.  The  stream  lines  of  the  combined  flow  for  the  northern  hemisphere 
at  the  500-mb  level. 
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